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Spin-orbit coupling (SOC), the interaction between the spin and momentum of a quantum particle, 
is crucial for many important condensed matter phenomena. The recent experimental realization of 
SOC in neutral bosonic cold atoms provides a new and ideal platform for investigating spin-orbit 
coupled quantum many-body physics. In this Letter, we derive a generic Gross-Pitaevskii equation 
as the starting point for the study of many-body dynamics in spin-orbit coupled Bose-Einstein 
condensates. We show that different laser setups for realizing the same SOC may lead to different 
mean field dynamics. Various ground state phases (stripe, phase separation, etc.) of the condensate 
are found in different parameter regions. A new oscillation period induced by the SOC, similar to 
the Zitterbewegung oscillation, is found in the center of mass motion of the condensate. 
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Spin-orbit coupling (SOC) for electrons plays a crucial 
role in many important condensed matter phenomena 
and applications, such as anomalous and spin Hall effects 
[H, topological insulator Q, spintronic devices etc. 
However, the observation of SOC physics in nature solid 
state systems are often hindered by the unavoidable dis- 
order and impurities. In this context, ultra-cold atomic 
gases provide an ideal platform for exploring novel SOC 
physics and device applications, owing to their unprece- 
dented level of control and precision in experiments. In 
ultra-cold atomic gases, SOC can be generated through 
the laser-atom interaction, which yields Abelian or non- 
Abelian gauge fields for atoms in the dressed state basis 
■ The recent broad interest in spin-orbit coupled cold 
atoms is mainly motivated by their remarkable applica- 
tions (3-[2fjj. ranging from the observation of Majorana 
fermions and the associated non- Abelian quantum statis- 
tics jL2jl(| , the design of atomtronics /spintronics devices 
B l9l l2l|. to the generation of magnetic monopoles Q , 
etc. An important difference between electrons and cold 
atoms is that, while electrons are fermions, ultra-cold 
atoms may be bosons, leading to novel spin-orbit physics 
that has not been explored in solids. 

A benchmark experiment along this research direction 
is the recent realization of the one-dimensional (ID) SOC 
for cold bosonic atoms [l5j] . which brings a completely 
new avenue for the study of the many-body dynamics 
of spin-orbit-coupled Bose-Einstein condensates (SOC- 
BECs). It is well-known that the starting point for the 
investigation of many-body dynamics of a BEC is the 
mean-field Gross-Pitaevskii (G-P) equation [221 ] . whose 
general formula is still lacked for the SOC-BEC. In this 
Letter, we derive a generic G-P equation for the SOC- 
BEC and investigate their ground states and collective 
excitations. Our main results are: 

(i) We find that the mean field dynamics depend on 
not only the SOC itself, but also the method to gener- 
ate it because the effective pseudospin states in different 
laser setups for the realization of the same SOC can be 
different superpositions of atom hyperfinc states with dif- 



ferent s-wave scattering lengths. This fact is now taken 
into account in the mean-field interaction terms in the 
pseudospin space in the G-P equation. We find that the 
mean-field interaction energy may depend not only on 
the density, but also the phase of the condensate (i.e., 
terms like + c.c), which, to the best of our knowl- 

edge, have not been explored in previous literature 0- 
|20( ■ We show that the phase dependent mean field terms 
can emerge for one type of laser setup (denoted as a com- 
plex system) , but vanish for another (denoted as a simple 
system), with both laser setups implementing the same 
Rashba SOC. 

(ii) We analyze the condensate wavefunctions for both 
systems in various parameter regions. In the strong inter- 
action and SOC region, there exist two distinct phases for 
the condensate density: Thomas-Fermi (TF) and stripe. 
While in the weak and medium interaction or SOC re- 
gion, spatial separation between two pseudospin compo- 
nents is observed. 

(iii) The low energy collective excitations in spin- 
orbit coupled BECs are investigated, for the first time, 
through the center-of-mass (COM) motion of the conden- 
sate when a sudden shift of the center of the harmonic 
trap is applied. We find a novel, shift direction and 
distance dependent, oscillation frequency in the COM 
motion induced by the SOC, similar as the Zitterbewe- 
gung (ZB) oscillation in the free space However, the 
new oscillation period is linearly proportional to the SOC 
strength (in contrast to the inverse proportion in the ZB 
oscillation [9]). It is also independent on the atom inter- 
action strength. Both oscillation period and amplitude 
are different for the simple and complex systems, even 
though they share the same physical parameters. 

Consider ultra-cold bosonic atoms confined in a quasi- 
two-dimensional (xy plane) harmonic trap with a tripod 
electronic level scheme (Fig. QJi). The atom dynam- 
ics along the z direction are frozen by a deep optical 
trap or lattice. The hyperfine ground states |1), |2), |3) 
are coupled with the excited state |0) using three blue- 
detuned lasers with the Rabi frequencies f2i, f^, and 
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Figure 1: (Color online) A tripod scheme for implementing 
Rashba spin-orbit coupling, (a) The atom-laser coupling, (b) 
The lowest three energy levels composed of hyperfine ground 
states |1), |2), |3). \Di) and ID2} are degenerate dark states, 
fi = ,710^1 + + I . (c,d) Two different laser configura- 
tions [a, for the tripod scheme. 

The Hamiltonian of the system can be written as 
H = J2i Hs( r i) + H int with the single particle Hamilto- 
nian 

H s (v)=p 2 /2m + V(r) + H I (r), (1) 

where V(r) — m (w^r 2 + w 2 z 2 ) /2 is the harmonic 
trapping potential with the corresponding trapping 
frequencies uj z and u±. Hj{r) = — HA |0) (0| — 
|0) (1| +0 2 |0) (2| +n 3 |0) (3\ +H.c.) describes the 
atom-laser interaction, and A is the detuning to the state 
|0). The diagonalization of the Hamiltonian Hj(r) yields 
two degenerate darks states \D\) — Yl a =i du{r)\a), 
\ D 2) = Ea=i /<*( r )l a ) and tw0 bright states (Fig. [lb), 
where the coefficients d a (r), f a (r) are determined by 
the laser parameters fli and A [23]. In experiments, a 
large detuning A is chosen to suppress the spontaneous 
emission of photons that heats the atom gas. The blue- 
detuned lasers are used to ensure the degenerate dark 
states are the ground states of the system to avoid the 
collision loss. In the subspace spanned by the pseudospin 
states |t) = l-Di) and \\) = \D 2 ), 

H s {r) = p 2 /2m + j(p x a y -p y a x ) + V(r), (2) 

where 7, the Rashba SOC strength, is equal to hkj2\f%m 
and (a/2 — l)hk/m for the two laser setups shown in 
Figs. QJ: and [TJi, respectively, k is the wavevector of 
the lasers. The s-wave scattering interaction between 
atoms can be written as H int = J2i<j E ai/ 8=i ffa^( r ? - 
r^) \a) - \(3)j (a\ i {j3\-, where rf is the position of the atom 
i in the hyperfine state \a) il g a p — Ai:h 2 a a p/m, a a/ 3 is the 
s-wave scattering length between atoms in the hyperfine 
states I a) and |/3). 

In the Hartree approximation, the many-body 
wavefunction of the bosonic system can be taken 
as the product of the single-particle wave function 
|*(ri,r 2 ,...,rjv)) = n£il*(r<)>, where |$(r)) = 
<&-|-(r) |Z?i (r)) + $^(r)|Z?2(r)) with the normalization con- 
dition Jdr(|$ t (r)| 2 + |<I^(r)| 2 ) = 1. The G-P equation 



can be obtained through the standard minimization of 
the mean-field energy functional E(&^, $4,) = (<I'|iJ| , I') 
with respect to $| and $| (22[. Because the scat- 
tering length is between atoms in different hyperfine 
states, the interaction energy (*P|.Hj„ t |*I>) should be eval- 
uated within the hyperfine state basis, i.e., |$(r)) — 
ELi^tW^W + $4.(r)/ Q (r)]|a>. With a straightfor- 
ward calculation for such energy functional minimization, 
we derive the G-P equation for a SOC-BEC 

ihd$/dt = if s $ + r$. (3) 

Here the two component wavefunction $ = ($-(-, §i) T in 
the pseudospin basis {|f) , The nonlinear term V = 

(li where r 1 =C 1 |$ t | 2 +C 2 |^| 2 +2Re(C 3 $^), 

r 2 = C 3 |$ T | 2 + C 4 |^| 2 + C 5 $* f $ h and T 3 = C 2 |$ T | 2 + 
C 6 |$_ l | 2 + 2Re(C 4 $|$ i ) C 3 , C 4 and C 5 terms, which 
are absent in previous study, originate from the linear 
superposition of hyperfine states for a pseudospin state. 
We see that not only the density, but also the relative 
phase between two components, play an important role 
on the dynamics of the BEC. Although our derivation for 
the mean-field interaction terms is based on the tripod 
schemes illustrated in Fig. [TJ it also applies to other 
schemes for gen erating Rashba SOC using more hyperfine 
states (e.g. [25j). 

For simplicity, we assume g a p = g\ for a ^ (3 and 
9a/3 = 9o for a — j3 with the corresponding scattering 
lengths ai and dQ. However the nonlinear G-P equation 
(|3"|) works for any g a p. Ci depend strongly on the laser 
configurations for implementing the Rashba SOC. Here 
we consider two laser setups for the tripod level scheme 
(Fig. [lji) that have been investigated previously in the 
literature [alia 23]. In the first setup (Fig. [TJ;), we have 
the SOC strength 7 = fc v /fi/12mo; ± , C* 3 = C 4 = C 5 = 0, 
Ci = C 6 = x(2a +4ai)/3, an d C 2 = x(4 ap + 2 ai )/3 (the 
simple system) with x — 2-/V yj2itm<jj z /h. Here we rescale 
the G-P equation ([3]) with the energy, time and length 
units huj±, wj 1 , and sj h/muj± respectively. The unit for 
7 is a velocity ^haj±/m. In the second setup (Fig. [TJl), 
7 = (V5 - l)ky/h/mu ± , d = C 6 = X [(12 - 8V2)a + 
{8V2 - 10)oi], C 2 = x[(24 - 16V2)a + (16^2 - 22)oi], 
C 3 = C 4 = -i X (7V2- 10)(ao - 01), and C 5 = x(4\/2 - 
6)(ao — di) (the complex system). 

We numerically solve the G-P equation ([3]) using the 
imaginary time evolution method and obtain the ground 
state of the condensate. In the simple system, there 
are two different types of phases in the region of strong 
SOC (7 > 1) and interaction (d > 1): the Thomas- 
Fermi (TF) phase (Fig. Hi) when C x > C 2 (equiva- 
lent to a\ > do) and the stripe phase when C± < C 2 
(Fig. (2t) [3. In the TF phase, the maximum densities 
of two components locate at the harmonic trap center, 
but the phase of the condensate varies like a plane wave 
exp (ik ■ r) (Fig. [3p) . In the stripe phase, the density for 
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Figure 2: (Color online) Condensate wavefunction $ in the 
simple system. (a,b) Density |$t| 2 ( a ) an d phase arg($-|-) (b) 
in the TF phase. Ci = 10, C2 = 6, 7 = 10. (c,d) In the stripe 
phase. Ci = 6, C 2 = 10, 7 = 10. (e,f) Density j$ t j 2 and 
|$;| 2 with Ci = 10, C 2 = 6, 7 = 1. (g,h) Density j$ t j 2 and 
|$J 2 with Ci = C 2 = 0,7 = 10. 



each pseudospin component forms a set of stripes (Fig. 

and two components are spatially separated. There 
is a sharp change of the condensate phase in and out- 
side the stripe region (Fig. [2}i) . The total density of two 
components has a TF distribution for both phases. The 
direction of the condensate phase variation V arg (<J> CT ) is 
spontaneously selected, along which the density distribu- 
tion is wider. We find that the density profile in each 
hyperfine state has the same structure (TF or stripe) as 
that in the pseudospin state, therefore the TF and stripe 
phases should be observable by directly measuring the 
density in each hyperfine state. 

In the medium SOC region (7 ~ 1), there exists spa- 
tial separation between the atom densities of two pseu- 
dospin components when C\ > C2 (Figs. [2]3, ET)- The 
separation can be understood from the spin-dependent 
force F = — 2 r y 2 (k x e z )a z generated by the Rashba 
SOC, where k is the momentum of atoms and along the 
condensate phase variation direction. F has the oppo- 
site directions for two pseudospins and is along the x 
direction when the phase variation is along the y direc- 
tion (Fig. Eb), leading to the spatial separation of two 
pseudospin components along the x direction (Figs. [2b, 
[2f). The separation due to the spin-dependent force is 
more transparent in the non- interacting region (Figs. [2tj, 
[2h) without involving the complexity from the interac- 
tion. However, the separation between two components 
decreases not only for very strong SOC (7^1, Fig. [5^), 
but also for very weak SOC (7 <C 1), which can be un- 
derstood based on F oc 7 2 and the zero separation for a 
regular spinor BEC without SOC. 

In the complex system, similar TF and stripe phases 
also exist in the strong SOC and interaction region. How- 
ever, significant differences between the simple and com- 
plex systems exist for the low energy collective excita- 
tions even though their condensate phases are similar. 
Here we consider the COM motion of the condensate 
subject to a sudden shift of the center of the harmonic 



trap. It is well-known that the COM motion of a BEC 
without SOC is a dipole oscillation whose frequency is 
the harmonic trap frequency and does not depend on the 
nonlinearity [26| . On the other hand, the interference 
between two Rashba spin-orbit energy bands yields the 
ZB oscillation [9] for a single particle in the free space 
with the oscillation period inversely proportional to 7 2 . 
In this Letter, we study the COM motion of a BEC in 
the presence of interaction, a harmonic trap, and Rashba 
SOC. 

The COM motion 

(r(i)) = y (|4> t (r,t)| 2 + |^(r,t)| 2 )rdr (4) 

of the condensate can be obtained by numerically solving 
the G-P equation © and is plotted in Figs. [3^ and [SJd 
(Figs. [3b, [3b!) for the TF phase in the simple (complex) 
system. In the simple system, when the shift of the cen- 
ter of the harmonic trap D is along the x direction (per- 
pendicular to the condensate phase variation V arg (3>f) 
direction y), the COM motion along the x direction is 
perfectly periodic with two periods (Fig. [3b,): one corre- 
sponds to the harmonic trap frequency, the other is much 
larger and linearly proportional to the SOC strength 7 
(Fig. 2]). The COM motion along the y direction is sim- 
ilar. However, when the shift of the harmonic trap is 
along the y direction, the COM motion along the y di- 
rection possesses only the harmonic trap frequency and 
the oscillation amplitude is strongly damped (Fig. 3b). 
The COM motion along the x direction vanishes. 

The different COM motions along different shifting di- 
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Figure 3: (Color online) The center of mass motion. Red 
lines: (a;); Blue lines: (y). (a,b) The simple system with the 
corresponding condensate wavefunction in Figs. [2k and [2b. 
The shift of the center of the harmonic trap D — 0.8. The 
shifts are along the x (a) and y (b) directions respectively. 
(c,d) The same COM motion as (a,b) but for the complex 
system. Parameters are the same as that in (a,b) except that 
new terms C3, C4, C5 are included in the condensate and 
COM motion caculations. 
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Figure 4: (Color online) Plot of the COM oscillation period 
T with respect to the SOC strength 7 for the simple system 
in the TF phase. Circles are from the numerical stimulation 
of the G-P equation ([3]). The lines are from Eq. ((5]). 

rections may be understood from the single atom dy- 
namics with the Rashba SOC. The spin-orbit coupled 
atoms have two bands with energies E± = h 2 k 2 /2m ± 
jk and corresponding wavefunctions <fi± = exp(ik • 
r)(l,±ie i ^ k ) T /v / 2, where tp k = arg(fc x + iky). The 
ground state of the atom 0_ stays at the potential min- 
imum located at kg = jm/h 2 with the direction of ko 
spontaneously selected. The shift of the harmonic trap 
corresponds to adding a momentum p into the ground 
state, leading to an initial state 4>i n i = exp(ip • r)</>_. 
The time evolution of the wavefunction can be written 
as ip{r,t) = J dk[A_ exp(iE-t)cj)- + A + exp(iE + t)<j>+], 
where A± = {<t>±\<j>i n i) ■ When p is perpendicular to the 
direction of ko (Fig. [2Jd), both A± are nonzero, and there 
is an interference between two spin-orbit coupled bands, 
leading to the SOC dependent oscillation of the COM 
motion (Fig. In contrast, when p is along the direc- 
tion of ko, one of A± must be zero and there is only one 
oscillation frequency due to the harmonic trap. 

We numerically calculate the new oscillation period T 
in Fig. [3^ for several different sets of interaction param- 
eters, and find that it does not depend on the interaction 
strengths C\ and C2 , indicating essentially single particle 
physics in this system. In Fig. |U we plot the dependence 
of T with respect to 7 for three harmonic trap shifts D. 
The numerical data can be well fitted with an analytic 
formula 

T = 2n{l + 1 /D). (5) 

We see that the period is linearly proportional to 7 and 
D^ 1 . In the strong SOC region (7 > 1), the strong SOC 
couples many harmonic oscillator states and tends to re- 
duce the single particle energy splitting, as confirmed by 
our numerical calculation, resulting in a large oscillation 
period. Note that this is very different from the ZB oscil- 
lation in the weak SOC region where the oscillation 
period is proportional to 7~ 2 because the energy split- 
ting between two spin-orbit coupled bands in this region 
is proportional to 7 2 , as confirmed in our numerical simu- 
lation. Finally, a smaller shift D leads to less excitations 
to the high energy states, yielding a larger COM oscilla- 
tion period. As 7 — > 0, T —> 2tt, the oscillation period for 
the harmonic trap frequency, as expected. In the weak 



and medium interaction or SOC region, our numerical 
results show that the oscillation amplitude decays with 
time due to the phase separation of the densities of two 
pseudospin components (Figs. [2^-[2h) that leads to the 
decoherence in the COM motion. 

In the stripe phase of the simple system, the COM mo- 
tion is always damped because of the decoherence orig- 
inating from the spatial separation of two spin compo- 
nents. Similar as the TF phase, the new oscillation period 
disappears (emerges) when the shift is along (perpendic- 
ular to) the condensate phase variation direction. 

The COM motion is strongly modified in the complex 
system. In Fig. [3fc and [3ji, we plot the COM motion in 
the TF phase. Clearly, the direction dependence of the 
new oscillation period is the same as that for the simple 
system. However, the oscillation period T is very differ- 
ent from that in the simple system (Fig. |3^) although the 
complex system has the same parameters (C±, C2, 7) as 
the simple one. A damping of the oscillation amplitude 
is also observed, in contrast to the perfect oscillation in 
the simple system. Furthermore, when the shift of the 
harmonic trap is along the condensate phase variation 
direction, the COM motion perpendicular to the shift di- 
rection emerges after a long time (Fig. [3ji), in contrast 
to the zero motion in the simple system (Fig. [3Jd). 

In experiments, we can choose the trapping frequency 
u±_ = 2tt x 20Hz, which yields the time unit wj 1 = 
8ms and length unit ah = y^h/mui^ = 2.4 /im. For Rb 
atoms, the spin-orbit coupling strength 7^1 — 10. The 
interaction strength C\, C2 can be tuned through the 
Feshbach resonance or by adjusting the atom number. 
With ao,ai ~ lOOas for Rb atoms, C\, C2 ~ 10 for 
N ~ 300, where 0,3 — 0.53 A is the Bohr radius. We 
consider a multiple layer system generated by an optical 
lattice along the z direction, therefore the total number 
of atoms is ~ 10 4 . We also confirm that the mean field 
dynamics are similar for larger interaction strength C\, 
C2 ~ 100, therefore the total number of atoms can be 
~ 10 5 . For an initial shift of the harmonic trap D = 2 
and 7 = 6, we find the center of mass oscillation period 
T ~ 200 ms, which is much shorter than the lifetime of 
the BEC and should be observable in experiments. 

In summary, we derive a generic G-P equation and in- 
vestigate the mean field dynamics for SOC-BECs. We 
emphasize that our G-P equation may serve as the start- 
ing point for the future study on the dynamics of SOC- 
BECs. Our predicted new oscillation period in the COM 
motion is observable in experiments and may provide 
a powerful tool for exploring the phase of the conden- 
sate wavefunction as well as the low energy excitations 
in SOC-BECs. 
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